#pragma once

#include <assert.h>
#include <QVector>
#include <QString>

#include "datainterface.h"

namespace PBVP {
	///
	class Reg3DVolumeData : public DataInterface {
	public:
		/// Default Constructor
		Reg3DVolumeData(unsigned int dim[3]);

		///
		Reg3DVolumeData(unsigned int dx, unsigned int dy, unsigned int dz);

		/// Copy Constructor
		Reg3DVolumeData(Reg3DVolumeData *pVolData);

		/// Destructor
		~Reg3DVolumeData();

		///
		DataSetType getDataSetType() const {
			return VOLUME_DATA;
		}

		unsigned int numOfVariables() const {
			return m_NumOfVar;
		}
		void setNumOfVariables(unsigned int n);

		///
		inline void getDimension(unsigned int dim[3]) const;
		inline unsigned int width() const { return m_dim[0]; }
		inline unsigned int height() const { return m_dim[1]; }
		inline unsigned int depth() const { return m_dim[2]; }
		inline void setWidth(unsigned int w) { m_dim[0] = w; }
		inline void setHeight(unsigned int h) { m_dim[1] = h; }
		inline void setDepth(unsigned int d) { m_dim[2] = d; }

		inline void getOrigin(float orig[3]) const;
		inline void setOrigin(const float orig[3]);

		inline void getSpan(float span[3]) const;
		inline float spanX() const { return m_span[0]; }
		inline float spanY() const { return m_span[1]; }
		inline float spanZ() const { return m_span[2]; }
		inline void setSpan(const float span[3]);
		inline void getMinMaxExt(float minExt[3], float maxExt[3]) const;
		inline void getMinMaxExt(double minExt[3], double maxExt[3]) const;
		inline void setMinMaxExt(float minExt[3], float maxExt[3]);
		inline void setMinMaxExt(double minExt[3], double maxExt[3]);
		inline double getVolume() const;

		DataType getDataType(unsigned int var = 0) {
			assert(var < numOfVariables());
			return m_types[var];
		}

		DataPtr getDataPtr(unsigned int var = 0) {
			assert(var < numOfVariables());
			return m_DataPtrs[var];
		}

		/// get the number of vertices
		inline unsigned int numOfVerts() const {
			return m_dim[0]*m_dim[1]*m_dim[2];
		}

		inline unsigned int numOfCells() const {
			return (m_dim[0]-1)*(m_dim[1]-1)*(m_dim[2]-1);
		}
		///
		inline double getValue(int nv, unsigned int var = 0) const;
		///
		inline double getValue(unsigned int i, unsigned int j, unsigned int k, unsigned int var = 0) const {
			//if (i >= 0 && i < m_dim[0] &&
			//	j >= 0 && j < m_dim[1] &&
			//	k >= 0 && k < m_dim[2])
			//{
				return getValue(index2vert(i, j, k), var);
			//}
		}

		/**
		* Get the function value at an arbitrary point
		* @note return 0 if the coordinates are outside of the bounding box
		*/
		double getValue(float coord[3], unsigned int var = 0) const;

		inline void setValue(double val, unsigned int nv, unsigned int var = 0);

		inline void setValue(double val, unsigned int i, unsigned int j, unsigned int k, unsigned int var = 0) {
			setValue(val, index2vert(i, j, k), var);
		}

		/// Get function values at eight corner of the cell (i, j, k)
		void getCellValues(double vals[8], unsigned int i, unsigned int j, unsigned int k, unsigned int var = 0) const;

		///
		void setDataArray(void *data, unsigned int var = 0, DataType type = FLOAT_TYPE);

		/**
		* Compute the gradient vector at all grid points.
		*/
		void calcGradient();

		/// Get the gradient at a grid point
		void getVertGrad(unsigned int i, unsigned int j, unsigned int k, double grad[3], unsigned int var = 0);

		/// Get the forward gradient at a grid point
		/**
		 *
		 * @note Forward derivative = (F(i+1)-F(i))/dx
		 */
		void getVertGradForward(unsigned int i, unsigned int j, unsigned int k, double grad[3], unsigned int var = 0);

		/// Get the backward gradient at a grid point
		/**
		 *
		 * @note Backward derivative = (F[i] - F[i-1])/dx
		 */
		void getVertGradBackward(unsigned int i, unsigned int j, unsigned int k, double grad[3], unsigned int var = 0);

		/// Compute upwind derivatives using quadratic ENO method
		/**
		 *
		 */
		void getQuadENODerivatives(unsigned int i, unsigned int j, unsigned int k, double plus[3], double minus[3], unsigned int var = 0);

		/// Get 2nd partial derivatives at a grid point
		/**
		 * @param par2 Returned partial derivatives. \
		 *		par[0] = Fxx, par[1] = Fxy, par[2] = Fxz, par[3] = Fyy, par[4] = Fyz, par[5] = Fzz
		 * The 2nd partial derivatives are computed using central difference
		 */
		void getVertParDeriv2(unsigned int i, unsigned int j, unsigned int k, double par2[6], unsigned int var = 0);

		void getCellGrads(unsigned int i, unsigned int j, unsigned int k, double grads[8][3], unsigned int var = 0);

		/// Get the gradient at a given point
		/**
		* @note return (0, 0, 0) if the coordinates are outside of the volume
		*/
		void getVertGrad(float x, float y, float z, double grad[3], unsigned int var = 0);

		bool hasGradient(unsigned int var = 0) const {
			return (m_Gradients[var] != NULL);
		}

		double getMinValue(unsigned int var = 0);
		double getMaxValue(unsigned int var = 0);

		/// Resample the data into an unsigned char array
		void resampleUChar(unsigned char ucArray[], unsigned int dim[3], unsigned int var = 0);

		/// Simple sampling by copying the data into a subvolume of the output array
		void resampleUCharSimple(unsigned char ucArray[], unsigned int dim[3], unsigned int var = 0);

		///
		unsigned int index2vert(unsigned int i, unsigned int j, unsigned int k) const {
			return (k*m_xy + j*m_dim[0] + i);
		}

		int index2vert(int idx3[3])
		{
			return(idx3[0] + idx3[1]*m_dim[0] + idx3[2]*m_dim[0]*m_dim[1]);
		}

		void vert2index(int id, int idx3[3])
		{
			idx3[0] = id % m_dim[0];
			idx3[1] = (id / m_dim[0]) % m_dim[1];
			idx3[2] = id / (m_dim[0]*m_dim[1]);
		}

		///
		unsigned int index2cell(unsigned int i, unsigned int j, unsigned int k) const {
			return( i + j * (m_dim[0]-1) + k *(m_dim[0]-1)*(m_dim[1]-1));
		}

		void cell2index(unsigned int id, unsigned int idx3[3])
		{
			idx3[0] = id % (m_dim[0]-1);
			idx3[1] = (id / (m_dim[0]-1)) % (m_dim[1]-1);
			idx3[2] = id / ((m_dim[0]-1)*(m_dim[1]-1));
		}

		/// overloaded function
		void cell2index(int id, int idx[3])
		{
			idx[0] = id % (m_dim[0]-1);
			idx[1] = (id / (m_dim[0]-1)) % (m_dim[1]-1);
			idx[2] = id / ((m_dim[0]-1)*(m_dim[1]-1));
		}

		/// Compute the coordinates of the ith vertex
		void getCoordinates(float coord[3], int id)
		{
			int idx[3];
			vert2index(id, idx);
			for(int i = 0; i < 3; i++) {
				coord[i] = m_orig[i] + idx[i]*m_span[i];
			}
		}

		void getCoordinates(float coord[3], int i, int j, int k)
		{
			coord[0] = m_orig[0] + i*m_span[0];
			coord[1] = m_orig[1] + j*m_span[1];
			coord[2] = m_orig[2] + k*m_span[2];
		}

		///
		void calcMinMax();

		/// Copy and append the data from another Reg3D structure
		/**
		 * The dimensions of the two Reg3D volumes must be same.
		 */
		void mergeData(Reg3DVolumeData *pVolData);

		/// Get the max absolute value among (i, j, k), (i+1, j, k), (i-1, j, k), (i, j-1, k), (i, j+1, k), (i, j, k-1) and (i, j, k+1)
		double getOneNeighborAbsMax(int i, int j, int k, unsigned int var = 0);

		/// Get the min absolute value among (i, j, k), (i+1, j, k), (i-1, j, k), (i, j-1, k), (i, j+1, k), (i, j, k-1) and (i, j, k+1)
		double getOneNeighborAbsMin(int i, int j, int k, unsigned int var = 0);

	protected:
		///
		void getSevenNeighorValues(unsigned int i, unsigned int j, unsigned int k, double vals[7], unsigned int dir /*= 0, 1, or 2 */, unsigned int var = 0);

		///
		void deriv2SevenValues(double vals[7], double d1[6], double d2[5], float dx);

		/// min and max values of the scalar function
		// double	m_fmin, m_fmax;
		/// dimension of data
		unsigned int m_dim[3];
		/// coordinates of the volume origin
		float	m_orig[3];
		/// span of each volume cell
		float	m_span[3];
		/// Number of variables to handle multi-variable volume data
		unsigned int  m_NumOfVar;

		QVector<DataType> m_types;
		QVector<DataPtr>  m_DataPtrs;
		QVector<double *>  m_Gradients;
		QVector<QString>  m_VarNames;

		// Min and Max of functional values
		QVector<double>	m_MinValues;
		QVector<double> m_MaxValues;

	private:

		unsigned int m_xy;
		static const int vertinfo[2][2][2];

		void init();

		void clear();

	};

	void Reg3DVolumeData::getDimension(unsigned int dim[]) const
	{
		dim[0] = m_dim[0];
		dim[1] = m_dim[1];
		dim[2] = m_dim[2];
	}

	void Reg3DVolumeData::getOrigin(float orig[3]) const {
		orig[0] = m_orig[0];
		orig[1] = m_orig[1];
		orig[2] = m_orig[2];
	}

	void Reg3DVolumeData::setOrigin(const float orig[3]) {
		m_orig[0] = orig[0];
		m_orig[1] = orig[1];
		m_orig[2] = orig[2];
	}

	void Reg3DVolumeData::getSpan(float span[3]) const {
		span[0] = m_span[0];
		span[1] = m_span[1];
		span[2] = m_span[2];
	}

	void Reg3DVolumeData::setSpan(const float span[3]) {
		m_span[0] = span[0];
		m_span[1] = span[1];
		m_span[2] = span[2];
	}

	void Reg3DVolumeData::getMinMaxExt(float minExt[3], float maxExt[3]) const
	{
		for (int i = 0; i < 3; i++)
		{
			minExt[i] = m_orig[i];
			maxExt[i] = m_orig[i] + (m_dim[i]-1)*m_span[i];
		}
	}

	void Reg3DVolumeData::getMinMaxExt(double minExt[3], double maxExt[3]) const
	{
		for (int i = 0; i < 3; i++)
		{
			minExt[i] = m_orig[i];
			maxExt[i] = m_orig[i] + (m_dim[i]-1)*m_span[i];
		}
	}

	void Reg3DVolumeData::setMinMaxExt(float minExt[3], float maxExt[3])
	{
		for (int i = 0; i < 3; i++)
		{
			m_orig[i] = minExt[i];
			m_span[i] = (maxExt[i] - minExt[i])/(m_dim[i]-1);
		}
	}

	void Reg3DVolumeData::setMinMaxExt(double minExt[3], double maxExt[3])
	{
		for (int i = 0; i < 3; i++)
		{
			m_orig[i] = (float)minExt[i];
			m_span[i] = (float)((maxExt[i] - minExt[i])/(m_dim[i]-1));
		}
	}

	double Reg3DVolumeData::getVolume() const {
		double len[3];
		for (int i = 0; i < 3; i++) {
			len[i] = m_span[i] * (m_dim[i]-1);
		}
		return (len[0]*len[1]*len[2]);
	}

	double Reg3DVolumeData::getValue(int nv, unsigned int var) const{
#ifdef _DEBUG
		assert(var < numOfVariables());
#endif
		DataType type = m_types[var];
		DataPtr data = m_DataPtrs[var];
		switch (type) {
		case UCHAR_TYPE:
			return data.uchar_ptr[nv];
		case USHORT_TYPE:
			return data.ushort_ptr[nv];
		case FLOAT_TYPE:
			return data.float_ptr[nv];
		case INT_TYPE:
			return (data.int_ptr[nv]);
		case DOUBLE_TYPE:
			return (data.double_ptr[nv]);
		}
		// should not come here
		assert(0);
		return 0;
	}

	void Reg3DVolumeData::setValue(double val, unsigned int nv, unsigned int var) {
		assert(var < numOfVariables());
		DataType type = m_types[var];
		DataPtr data = m_DataPtrs[var];
		switch (type) {
		case UCHAR_TYPE:
			data.uchar_ptr[nv] = (unsigned char)val;
			break;
		case USHORT_TYPE:
			data.ushort_ptr[nv] = (unsigned short)val;
			break;
		case FLOAT_TYPE:
			data.float_ptr[nv] = val;
			break;
		case INT_TYPE:
			data.int_ptr[nv] = (int)val;
			break;
		case DOUBLE_TYPE:
			data.double_ptr[nv] = val;
			break;
		}
	}

}
